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model  is  suitable  for  incorporation  into  finite  difference  and  finite  element  stress  wave 
codes.  Note  that  numerical  simplicity,  efficiency,  and  stability  are  essential  characteristics 
of  a  model  which  is  to  be  of  practical  utility. 

This  research  has  also  been  supported  by  the  Department  of  Energy,  under  LLNL 
Contract  No.  B 157346,  and  LLNL  UC  Master  Task  Agreement  Task  #4. 


DTK  QUALITY  INSPECTED  8 


Accesion  For 

NTIS  CRA&I 
OTIC  TAB  0 

Unannounced  □ 

Justification 


By _ 

Distribution  / 


Availability  Codes 


Dist 

Avail  and/or 

Special 

11 
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1.  OBJECTIVES 

Propagation  of  seismic  waves  in  the  nearfield  where  rock  rheology  is  demonstrably 
nonlinear  raises  unique  difficulties.  Nonlinearity  arises  primarily  in  two  forms  at 
intermediate  to  large  strains:  (1)  nonlinear  elasticity,  and  (2)  amplitude-dependent 
attenuation,  which  is  a  well  documented  behavior  at  intermediate  strains  and  low 
confining  pressure.  The  proper  representation  of  nonlinear  constitutive  equations  for  rocks 
in  this  regime  is  a  potentially  important  ingredient  of  quantitative  source  models. 

Stress  wave  propagation  and  attenuation,  in  rocks  and  soils,  show  evidence  of 
significant  nonlinearity  at  strain  amplitudes  as  low  as  10~6,  leading  in  particular  to  an 
amplitude  dependence  of  the  apparent  Q,  most  likely  associated  with  friction  along 
microcracks  and  joints  [e.g.  Boitnott,  1992].  Recent  quasistatic  laboratory  testing  of  rock 
at  low  strain  has  permitted  detailed  high-quality  observations  of  cusped  hysteresis  loops 
in  this  regime.  These  issues  have  been  recently  reviewed  by  Minster  et  al.  [1991]  and 
summarized  by  Martin  and  Minster  [1992].  Nonlinear  wave  propagation  in  geological 
materials  has  also  been  observed  and  modeled  in  a  different  context  by  Bonner  and 
Wannamaker  [1991],  and  by  Johnson  et  al.  [1991].  Our  objective  is  to  identify  and 
validate  a  rheological  model  (constitutive  equation)  for  rocks,  valid  at  moderate  strains, 
that  explains  satisfactorily  these  various  observations,  and  is  appropriate  for  incorporation 
in  numerical  source  and  wave  propagation  codes,  and  apply  the  rheological  model  to 
improve  our  understanding  of  seismic  source  physics. 

We  have  shown  in  previous  work  that  nonlinear  one-dimensional  wave  propagation  can 
result  in  spectral  distortions  at  all  wavelengths,  but  that  this  effect  is  strongly  pulse-shape 
dependent,  and  therefore  call  for  a  3-D  capability  [Minster  et  al,  1991].  More  recently, 
we  have  found  that  our  use  of  an  approximate  description  of  the  phenomenological*' 
behavior  of  rocks  in  the  nonlinear  regime  is  flawed  insofar  as  it  is  not  able  to  simulate  new 
high-quality  laboratory  observations  of  hysteresis  loops  in  both  Sierra  White  granite  and 
Berea  sandstone  [Day  et  al.,  1992].  Ultimately,  a  complete  description  of  the  rheology  in 
terms  of  a  thermodynamically  valid  constitutive  equation  is  really  what  should  be  used  in 
numerical  simulations,  if  it  can  be  developed  and  validated  experimentally. 

2.  Nonlinear  Wave  propagation  and  Attenuation 

Our  earlier  numerical  modeling  of  nonlinear  attenuation  in  the  intermediate  strain 
regime  used  viscoelastic  theory  as  its  point  of  departure  [Minster  and  Day,  1986;  Minster 
et  al.  1991].  We  review  that  approach  here,  in  order  to  highlight  its  analogies  as  well  as 
its  contrasts  with  the  new  approach  proposed  in  Section  4.  A  further  reason  for  reviewing 
the  numerical  approach  to  viscoelasticity  is  that,  to  be  acceptable,  a  model  for  the 
nonlinear  intermediate  strain  regime  should  be  well  behaved  in  the  low  strain  limit. 
Thus,  it  will  be  desirable  to  develop  a  numerical  wave  propagation  treatment  which 
reduces  to  linear  viscoelasticity  in  the  small  amplitude  limit. 

2.1  Broad-band  approximation  of  amplitude-dependent  attenuating  rheology 

Based  on  a  suite  of  one-dimensional  simulations  of  nonlinear  wave  propagation 
problems  Minster  et  al.  [1991]  concluded  that  a  simple  model  in  which  QmI  is  simply 
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assumed  to  be  proportional  to  strain  amplitude  can  explain  the  shape  distortion  of  Lorentz 
peaks  observed  in  the  laboratory  at  moderate  strains,  and  the  apparent  superposability  of 
simple  pulses  even  in  the  nonlinear  regime.  They  also  concluded  that,  in  contrast  to  linear 
Q  models  for  which  the  spectrum  of  the  "Q  operator”  tends  to  unity  at  low  frequencies,  a 
nonlinear  rheology  may  lead  to  significant  spectral  distortions  at  all  frequencies,  and 
energy  losses  can  be  substantial  even  «t  wavelengths  long  compared  to  the  propagation 
distance.  Thus,  even  though  nonlinear  rheology  is  only  relevant  within  a  limited  distance 
from  a  seismic  source,  this  raises  the  possibility  that  the  far  field  source  spectrum  can  be 
affected  to  some  degree  at  all  frequencies,  including  those  pertinent  to  regional  phases  and 
teleseismic  body  waves. 

Those  results  were  based  on  an  attenuation  model  described  by 

Qr^Qd  +  re  ,  (2.i) 

where  e  is  the  strain  amplitude,  y  is  a  material  constant,  and  Q'J  represents  a  linear 
anelastic  term  controlled  by  mechanisms  that  mask  the  nonlinear  ones  at  low  strain.  This 
form  of  amplitude-dependence  describes  well  the  bulk  of  laboratory  evidence  accumulated 
to  date.  Nonlinear  wave  propagation  simulations  were  conducted  in  two  steps.  First,  we 
used  the  Padd  approximant  method  of  Day  and  Minster  [1984]  to  convert  the  stress-strain 
relation  of  a  linear,  anelastic  solid,  with  frequency-independent  Q,  into  differential  form. 
An  absorption  band,  with  Q  nearly  constant  at  Qo,  and  with  minimum  and  maximum 
relaxation  times  X\  and  Xi,  respectively,  yields  the  following  relation  between  stress 

history,  o(t)  and  strain  history,  e(t) 


deit')  , 


where  Mu  is  the  unrelaxed  modulus.  We  showed  that  (2.2)  can  be  approximated  by 


(2.2) 


<Xt)  =  Mu  |  £(/)-£  £(f) 

i=i 

where  the  6 's  are  relaxation  terms  governed  by  the  n  linear  equations 


(2.3) 


(2.4) 


The  constants  v,  and  w,  which  depend  on  the  order  of  approximation,  n,  are  given  by 
Day  and  Minster  [1984],  who  also  show  that  the  operator  defined  by  (2.3)  and  (2.4) 
converges  to  the  exact  result  (2.2)  as  n  increases.  The  second  step  is  to  generalize  (2.4) 
by  introducing  a  linear  dependence  of  Qo  on  strain  amplitude  according  to  (2.1): 

+  v,£  =  T]  w{Qf  +r|e(f)i)  £(/)  .  (2.5) 

Then,  (2.3)  and  (2.5)  constitute  the  stress-strain  equations  for  our  one-dimensional  finite 
difference  simulations. 

All  differential  operators  generated  by  this  procedure  can  be  guaranteed  to  be  causal, 
stable,  and  dissipative.  However,  that  the  method  performs  rather  poorly  when  the 
absorption  band  is  much  broader  than  the  calculational  pass  band,  that  is,  the  interval 
between  the  maximum  and  minimum  frequencies  resolvable  by  the  numerical  method. 
For  example,  the  finite  difference  method  is  limited  to  the  frequency  band  from  l/ndt  to 
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roughly  1/m#,  where  St  is  the  time  step,  n  is  the  total  number  of  time  steps  computed,  and 
m  is  the  number  of  time  steps  associated  with  the  minimum  resolvable  wavelength;  m  is 
typically  of  the  order  of  20,  and  n  may  be  up  to  several  thousand  for  large  two- 
dimensionai  calculations.  We  have  devised  a  simple  extension  of  the  method  which 
renders  it  suitable  for  broad  absorption  bands,  without  compromising  its  analytical  and 
numerical  simplicity.  Using  the  Laplace  transform  in  ^-multiplied  form,  we  reduce  the 
stress-strain  relation  to  its  operational  form: 

o(s)  =  M{s)e(s)  .... 

_  * 

Note  that  the  operational  modulus  M  has  the  same  dimensions  as  the  step  response  M. 
The  unrelaxed  modulus  Mu,  the  relaxed  modulus  Mr,  the  modulus  defect,  SM,  and  the 
normalized  relaxation  function  <j>,  are  given  by 


Mu  =  M(0)  =  M(oc)  ^ 

(2.7) 

Mr  =  A/(oo)  =  M(0) 

* 

(2.8) 

8M  =  MU-  Mr 

(2.9) 

M{t)  =  Mr  +SM  <p{t ) 

(2.10) 

We  represent  the  relaxation  function  in  terms  of  a  relaxation  spectrum  0, 

(P(t)  =  j  0(/n  t)  expi-t/t )  diln  t) 

(2-11) 

resulting  in  the  following  integral  expression  for  the  operational  modulus: 


M(s)  =  MU-  5M 


/s 

(pjp)dp 
s  +p 


(2.12) 


where  <p(p)  =  <2>(/n  T  '*)  .  We  may  now  partition  of  the  p  integral  into  3  regimes, 
separated  by  low-frequency  cutoff  pmin  and  high-frequency  cutoff  pmax) 


M(s)  =  Mu  -  SM  (/j  +I2  +  I3), 


[ Pmin 

r  ~ 

1  <kp)dp 

(Pip)  dP 

<Pip)dp 

1  s+p 

2  1  s  +  P 

3  1  s+p 

1 0 

J  Pmin 

J  Pmat 

(2.13) 

(2.14,15,16) 


The  interval  (pmtn,  pmax)  is  prescribed  to  coincide  with  the  calculational  pass  band.  The 
middle  partition,  I2,  is  replaced^ith  an  n^  order  Pad6  approximant,  as  before,  but  with  the 
support  interval  of  interval  of  0,  (t^1.  Tj1),  replaced  by  the  interval  (pmjn,  Pmax)  Then  lj 
and  I 3  are  approximated  by  Taylor  series  about  °o  and  0,  respectively.  Laplace  inversion 
leads  to  a  representation  of  total  stress  as  a  sum  of  n+2  internal  variables,  each  of  which 
satisfies  a  first  order  differential  equation.  This  permits  us  to  model  broad  absorption 
bands  efficiently  and  with  much  better  accuracy  than  before. 


3.  Laboratory  Observations  of  Hysteresis  loops 

In  order  to  validate  the  nonlinear  models  described  above,  we  have  conducted 
simulations  of  hysteresis  loops  measured  at  several  strain  amplitudes  in  uniaxial  tests  on 
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Sierra  White  granite  and  Berea  sandstone.  These  data  have  been  collected  by  New 
England  Research  Inc.,  and  have  been  kindly  made  available  to  us  by  Drs.  R.  Martin,  R. 
Haupt,  and  G.  Boitnott.  As  reported  by  Day  et  al.  [1992],  these  simulations  brought  to 
light  a  serious  shortcoming  of  our  approach,  namely  that  it  does  not  produce  the  correct 
loop  shapes  when  the  strain  amplitude  is  increased  into  the  nonlinear  regime.  Various 
modifications  of  our  general  approach-all  resulted  in  failure,  pointing  to  the  need  for  a 
completely  different  treatment  of  the  rheology,  dealing  intrinsically  with  the  nonlinearity. 

Laboratory  stress-strain  curves  under  cyclic  loading  characteristically  exhibit  the 
following  features  which  a  successful  model  must  emulate: 

•  Hysteresis  occurs,  implying  energy  loss,  and  the  effective  Q  characterizing  this 

dissipation  is  strain-amplitude  dependent. 

•  The  hysteresis  loops  are  cusped  at  reversal  points,  rather  than  elliptical  (as  would 

typify  linear  anelastic  behavior). 

•  No  yield  surface  is  evident  in  the  loading  curves,  at  least  for  strains  up  to  about  10-4. 

•  Upon  reversal  of  strain  path,  the  tangent  modulus  is  roughly  equal  to  the 

instantaneous  elastic  modulus. 

Typical  raw  laboratory  data  in  the  form  of  stress  and  strain  histories  are  often  rather 
noisy,  and  require  filtering.  Simple  low-pass  filtering  smoothes  the  cusps,  thereby 
masking  the  onset  of  nonlinear  behavior,  and  affecting  the  measurement  of  moduli  at  and 
near  the  cusps.  We  have  therefore  developed  a  technique  to  filter  separately  the  loading 
and  unloading  portions  of  the  loops.  It  relies  on  the  construction  of  a  longer  time  series 
out  of  a  half-loop — that  is  a  portion  of  stress-strain  history  between  two  reversals,  in 
which  both  stress  and  strain  are  monotonic.  This  is  done  by  extending  it  in  both 
directions  with  versions  of  itself,  rotated  by  ±n  about  its  end  points.  The  extended  time 
series  is  then  de-meaned,  de-trended,  and  low-pass  filtered  using  a  phaseless  filter  to 
avoid  introduction  of  a  phase  shift.  The  filtered  version  is  truncated  to  the  original  length 
after  restoring  trend  and  mean,  and  the  stress-strain  path  reconstructed  by  concatenation 
of  filtered  segments.  The  rotation  by  ±n  of  the  extensions  has  the  advantage  of 
preserving  the  continuity  of  the  time  series  and  its  derivative,  thereby  limiting  undesirable 
end  effects.  It  is  important  to  avoid  introducing  cusps  into  hysteresis  loops  when  they  are 
not  present,  and  to  avoid  smoothing  through  a  cusp  or  changing  its  angle,  when  one  is 
present.  Too  strong  a  filter  will  change  the  slope  near  the  end  points,  and  therefore 
introduce  a  fictitious  cusp.  Misapplication  of  the  technique  is  detectable  because  this  will 
create  overlaps  or  gaps  between  successive  segments.  The  technique  gives  very 
satisfactory  results  for  noise  levels  as  large  as  10  percent  as  we  have  been  able  to  verify 
using  synthetic  loops  contaminated  with  additive  noise.  This  approach  facilitates 
considerably  the  estimation  of  the  tangent  modulus,  particularly  near  the  loop  ends  where 
noise  contamination  is  most  worrisome. 

Several  of  the  features  described  above  are  evident  in  Figures  1  and  2,  which  show 
selected  sequences  of  hysteresis  loops  in  Berea  Sandstone  and  in  Sierra  White, 
respectively,  under  uniaxial  stress.  In  both  instances,  the  strain  and  stress  time  histories 
are  shown  in  the  top  frame,  followed  by  the  corresponding  stress-strain  paths  (hysteresis 
loops)  on  an  expanded  scale,  after  removal  of  the  mean  slope,  in  order  to  emphasise  the 
key  nonlinear  characteristics;  the  bottom  frame  illustrates  the  dependence  of  the  tangent 
modulus  on  strain.  In  the  filtered  data,  the  cusped  nature  of  the  reversal  points  is  evident. 
Also  evident  is  the  near-equality  of  the  initial  loading  and  unloading  slopes.  The  results 
illustrate  clearly  the  non-elliptical  (nonlinear)  character  of  the  hysteresis  loops  at  such 
moderate  strain  levels,  and  also  bring  out  clearly  the  strain  hardening  which  causes  the 
loops  to  show  upward  concavity. 


4 


Strain  (miooanin)  number*  indicate  Kpnents 

Figure  1:  (a)  Strain  (solid  line)  and  stress  (dotted  line)  time  histories  for  a  uniaxial  stress 
test  on  Sierra  White  granite,  (b)  corresponding  stress  and  strain  path,  after  smoothing; 
raw  data  are  indicated  by  the  dots,  and  the  mean  slope(modulus)  of  the  hysteresis  loops 
has  been  removed  to  emphasize  the  nonlinear  characteristics,  (c)  Young’s  modulus 
dependence  on  strain  for  this  test.  Numbers  indicate  the  various  segments  in  the  loading 
history,  separated  by  path  reversals. 
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Strain  (microwain)  numbers  indicate  segments 

Figure  2:  (a)  Strain  (solid  line)  and  stress  (dotted  line)  time  histories  for  a  uniaxial  stress 
test  on  Berea  Sandstone,  (b)  corresponding  stress  and  strain  path,  after  smoothing;  raw 
data  are  indicated  by  the  dots,  and  the  mean  slope(modulus)  of  the  hysteresis  loops  has 
been  removed  to  emphasize  the  nonlinear  characteristics,  (c)  Young’s  modulus 
dependence  on  strain  for  this  test.  Numbers  indicate  the  various  segments  in  the  loading 
history,  separated  by  path  reversals. 


4.  Endochronic  Constitutive  Model 


4.1  Introduction 

A  successful  rheological  model  should  be  capable  of  matching  the  features  seen  in 
Figures  1  and  2,  which  are  often  difficult  to  see  in  the  raw  data,  but  it  should,  as  much  as 
possible,  avoid  introducing  a  large  number  of  additional  model  parameters  for  this 
purpose.  The  endochronic  model  described  in  this  paper  offers,  we  think,  a  very 
promising  solution.  In  light  of  the  observations  outlined  in  the  previous  section,  we  have 
adopted  an  approach  to  modeling  the  moderate  strain  regime  which  departs  sharply  from 
viscoelastic  models,  yet  retains  much  of  the  computational  simplicity  described  in 
Section  2. 

4.2  Background  to  the  Endochronic  Formulation 

From  the  outset,  we  consider  the  class  of  constitutive  models  known  as  simple 
materials.  With  this  restriction,  the  stress  at  a  point  depends  only  on  the  strain  history  at 
that  point  (not,  for  example,  on  strain  gradients),  i.  e.. 


o(t)  =  F[e{t  ),0<f  <t]  (41) 

where  F  is  a  functional  relating  the  stress  a  to  the  strain  history  e{t).  For  example:  if 
F  is  linear  and  time  invariant,  equation  (1)  reduces  to  a  convolution,  and  we  have 
the  usual  formulation  of  viscoelasticity: 

0(t)  =  M(t)  *  e(t)  '  (4  2) 

This  restriction  combined  with  rate  independence  constitute  sufficient  conditions  to 
ensure  preservation  of  cube  root  scaling.  To  specialize  4.1  for  a  rate-independent  simple 
material,  we  express  the  strain  history  in  terms  of  the  strain  path  length  £.  Then 


0(0  =  F[e(D4.0<£(r  ')<««] , 

The  concept  of  rate.-independence  implies  that  there  is  no  dependence  of  the 
rheology  on  the  rate  £  : 


(4.3) 


o(0  =  F[c(|),  0<S(/')<£(0] 


where 


(4.4) 


(4.5) 


d t;  =  (de:g:de)m 

In  other  words,  £  is  the  strain  path  length,  measured  in  terms  of  the  metric  g. 

4.3  The  Endochronic  Material  Model 

Following  Valanis  and  Read  [1979],  we  consider  me  special  case  in  which  F  is  linear 
and  shift-invariant  in  the  plastic  strain  path  length  z 


where 


is  the  plastic  strain  increment.  The  linear,  shift-invariant  assumption  guarantees  that  we 
can  write  a  as  a  convolution  over  z:  that  is: 


dz  =  0 d&.g:d6)m 

(4.6) 

i 

■X3 

II 

(4.7) 

2  \L 

G(t)  =  K(z )  * 

dz 


(4.8) 


If  the  kernel  K{z)  is  chosen  to  have  an  integrable  singularity  at  z  =  0,  then  all  the  features 
noted  above  are  realized: 


K(z)  ocZa>0<CC<l'  (4  9) 

It  is  the  singular  behavior  of  the  kernel  that  insures  that  loading  and  unloading  at  reversal 
points  occurs  with  stress-strain  slope  equal  to  the  elastic  modulus.  Furthermore,  as 
demonstrated  below,  we  have  been  able  to  show  that  the  singular  kernel  ensures  power 
law  dependence  of  Q' 1  on  strain  amplitude,  in  accordance  with  experimental 
observations  cited  previously. 


4.4  Amplitude-Dependence  of  Q  in  the  Endochronic  Model 

The  power  law  amplitude  dependence  of  Q'  ^  is  derived  by  noting  that 


Restricting  treatment  to  uniaxial  loading, 


(4.10) 


so,  in  terms  of  the  maximum  plastic  strain  zm. 


=  1 


(4.11) 


Ooc  ( Zm)  . 

From  the  definition  of  Q~  *  in  terms  of  the  area  of  the  hysteresis  loop,  we  obtain: 

Q- 1  oc  amZm 


and  from  4. 13  and  4. 14,  we  obtain 


Qa  ~  . 


(4.12) 

(4.13) 

(4.14) 


Note  that  for  a  =  1/2,  we  have  an  approximately  linear  dependence  on  strain  amplitude, 
in  agreement  with  a  large  body  of  laboratory  observations.  The  endochronic  model  thus 
appears  capable  of  emulating  laboratory  observations  of  hysteretic  behavior,  as  well  as 
amplitude  dependence  of  attenuation  at  moderate  strain  amplitudes. 


4.5  Computational  approach 

To  be  useful  for  numerical  simulations,  the  convolution  form  4.8  must  be  converted  to 
a  differential  constitutive  equation.  Since  the  endochronic  model  has  a  formal  structure 
similar  to  linear  viscoelasticity,  we  can  carry  out  this  conversion  in  a  manner  analogous 
to  that  used  in  Section  2.  That  is,  we  first  Laplace  transform  4.8, 


Ot»  =  sK(s)6(s)  (415) 

We  approximate  K(s)  by  a  rational  function  Kn(s),  where  n  is  the  order  of  the 
denominator.  Then,  we  develop  Kn(s )  as  a  partial  fraction  expansion 


*"(*)  =  2 

£S+VJ' 


(4.16) 


where  Vj  and  A j  are  the  poles  and  residues,  respectively,  of  Kn(s). 

Transformation  back  to  the  z  domain  yields  the  following  system  of  differential  equations 
for  a: 
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(4.17) 


Ofc)  =  £Cj 


Equations  4.6,  4.7,  and  4.17  form  the  set  of  constitutive  equations  which  are  solved 
numerically.  A  convenient  numerical  scheme  for  the  solution  of  such  a  system  is  given 
by  Murakami  and  Read  [1989],  a.:d  we  have  successfully  implemented  that  scheme  to 
compute  the  numerical  results  shown  in  the  next  section. 


5.  Application  to  Laboratory  and  field  Data 

To  validate  the  use  of  the  endochronic  model,  we  have  conducted  simulations  of 
hysteresis  loops  measured  at  several  strain  amplitudes  in  uniaxial  tests  on  Sierra  White 
granite  and  Berea  sandstone  data  collected  by  New  England  Research  Inc. 

Figure  3  shows  a  comparison  of  such  simulations  with  three  observed  loops  in  Sierra 
White,  at  stress  levels  of  3,  6,  and  12  bars,  respectively.  The  numerical  simulations 
(which  are  symmetrical  in  stress  and  strain  histories)  match  well  the  overall  character  of 
the  observations;  this  includes  in  particular  the  increase  in  attenuation  with  increasing 
strain. 

Figure  4  shows  a  similar  comparison  for  Berea  sandstone,  for  which  the  attenuation 
levels  are  much  higher,  as  illustrated  by  the  loop  areas.  With  respect  to  attenuation  and 
loop  shape,  the  comparison  is  quite  favorable.  The  theoretical  loops  simulate  the 
amplitude  dependence  of  Q  and  the  non-elliptical  loop  shapes,  including  cusps  at  the 
ends.  The  mean  slope  decreases  somewhat  more  rapidly  with  increasing  strain  amplitude 
for  the  experimental  loops  than  it  does  for  the  theoretical  loops.  Further  improvement  of 
the  model  fit  to  the  Berea  Sandstone  data  can  be  obtained  by  introducing  variations  into 
the  shape  of  the  kernel  function. 

We  have  also  calculated  Q'*  as  a  function  of  strain  amplitude  for  both  the  Sierra  White 
Granite  and  Berea  Sandstone  models.  The  Berea  model  produces  Q"1  nearly  proportional 
to  strain  amplitude  over  the  full  range  shown.  The  present  Sierra  White  model  also 
exhibits  a  strong  amplitude  dependence  of  {T*.  although  it  departs  slightly  from  the 
expected  linear  dependence  of  Q~ 1  on  strain  amplitude.  This  is  probably  as  a  result  of 
approximations  introduced  in  our  current  expansion  of  the  singular  kernel  function. 

The  singrlar  kernel  endochronic  model  reproduces  several  key  nonlinear  phenomena  ^ 
associated  with  rock  hysteresis  at  moderate  strain.  The  approach  represents  a  substantial 
improvement  over  earlier  attempts  to  simulate  amplitude-dependent  attenuation  using 
variants  of  viscoelasticity.  Although  purely  phenomenological,  the  endochronic  approach 
has  the  decided  advantage  that  it  readily  reduces  to  a  set  of  relatively  simple  differential 
equations  which  are  easily  solved  numerically.  All  numerical  results  reported  here  were 
obtained  by  solving  this  system  of  differential  equations  numerically.  Exactly  the  same 
algorithm  can  be  applied  to  compute  stress-strain  behavior  in  numerical  wave 
propagation  codes. 


9 


Sierra  WVta  Granite  SWTIO  dao(_J  rywhrae*-) 


Figure  3:  Comparison  of  endochronic  simulations  with  three  observed  loops  in  Sierra 
White,  at  stress  levels  of  3,  6,  and  12  bars,  respectively.  The  numerical  simulations 
(which  are  symmetrical  in  stress  and  strain  histories)  match  well  the  overall  character  of 
the  observations;  this  includes  in  particular  the  increase  in  attenuation  with  increasing 

strain.  Additional  simulations  for  larger  stress  amplitudes  show  that  Q*1  continues  to 
increase  at  large  strain  levels. 
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Figure  4:  Similar  comparison  (see  Figure  3)  for  Berea  sandstone,  for  which  the 
attenuation  levels  are  much  higher,  as  illustrated  by  the  loop  areas.  Again,  the 
comparison  is  quite  favorable,  including  the  amplitude  dependence  of  Q,  and  the  non 
elliptical  loop  shapes,  with  apparent  cusps  at  the  ends.  It  should  be  emphasized  that, 
unlike  many  nonlinear  models,  the  model  used  in  these  simulations  depends  only  on  a 
small  number  of  parameters,  once  the  kernel  singularity  is  specified. 
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As  stated  before,  a  1-D  nonlinear  numerical  model  is  not  easily  generalized  to  3-D. 
For  stress  wave  propagation,  there  is  much  more  to  it  than  merely  including  geometrical 
spreading,  because  the  rheology  itself  is  amplitude-dependent.  Applications  of  this  class 
of  algorithms  to  the  interpretation  of  seismological  data  collected  in  the  field  require 
therefore  development  and  validation  of  a  full  3-D  wave  propagation  capability.  Intuition 
fails  us,  or  is  even  misleading  for  nonlinear  situations,  so  that  it  makes  little  sense  to 
develop  such  a  capability  until  the  model  has  been  fully  verified  on  laboratory  data  in  the 
1-D  situation.  We  will  therefore  defer  such  effort  until  later. 

6.  Conclusions 

Existing  theory  for  seismic  wave  propagation  is  almost  exclusively  linear,  despite 
abundant  experimental  and  observational  evidence  of  nonlinear  phenomena  in  earth 
materials  for  strain  levels  exceeding  about  10'6.  Our  long-term  goal  is  to  make  a 
significant  contribution  to  filling  this  gap  in  seismic  theory.  Our  approach  promises  to 
provide  a  stable  and  efficient  algorithm  for  numerical  simulation  of  nonlinear  wave 
propagation  at  intermediate  strain  levels,  with  computational  requirements  only  modestly 
exceeding  those  of  linear  viscoelasticity.  The  method  should  enable  numerical  modeling 
to  better  account  for  near-source  nonlinear  phenomena,  which  in  turn  will  improve  our 
understanding  of  source  physics  for  both  earthquakes  and  buried  explosions. 

In  particular,  the  importance  of  nonlinearity  in  the  intermediate  strain  regime  for 
detection,  identification,  and  yield  estimation  of  underground  explosions  remains  a 
significant  unresolved  issue.  Specific  phenomena  which  should  be  considered  include  ( 1 ) 
the  effects  of  nonlinear  attenuation  on  surface  reflections  (e.g.,  “depth”  phases  such  as  pP 
and  pS),  both  for  sources  sufficiently  shallow,  so  that  the  nonlinear  regime  extends  to  the 
free  surface,  and  for  the  case  of  strongly  attenuating  surface  layers  (soils);  (2)  effects  of 
nonlinear  attenuation  on  the  efficiency  of  high  frequency  cavity  decoupling;  and  (3)  the 
effect  of  nonlinear  attenuation  on  the  spectral  characteristics  of  regional  seismic 
recordings  from  both  shallow  and  overburied  explosions.  For  example,  Taylor  and 
Randall  [1989]  have  identified  systematic  spectral  differences  between  regional 
seismograms  from  shallow  explosions  and  overburied  explosions  at  NTS.  The  spectra  of 
regional  phases  play  an  important  role  in  event  identification  in  the  context  of  a 
nonproliferation  or  even  a  CTBT  treaty,  and  it  is  thus  very  important  to  establish  the 
physical  origin  of  such  spectral  differences.  Our  work  is  aimed  ultimately  at 
understanding  such  near-source  effects.  •> 

This  project  is  being  undertaken  in  coordination  with  experimental  work  at  New 
England  Research  Corporation  (NER).  Our  modeling  results  will  be  used  to  help  guide 
the  design  of  subsequent  NER  experiments.  Those  model  parameter  which  are  found  to 
be  most  critical  in  controlling  the  seismic  signature  of  explosions  should  be  identified  and 
targeted  for  experimental  study.  Such  collaboration  has  already  been  initiated,  and  all  the 
data  sets  shown  in  this  report  have  been  made  available  as  a  result  of  it.  Experiments 
conducted  by  NER  to  date  have  focused  on  uniaxial  stress  geometries.  However,  shear 
attenuation  is  most  important  in  the  Earth,  so  that  future  experiments  in  torsion  are  of 
particular  importance,  as  well  as  experiments  highlighting  the  effects  of  pore  fluids  and 
saturation,  which  are  essential  at  low  strains.  We  are  particularly  concerned  with  the 
ability  of  the  endochronic  model  to  accommodate  such  effects,  in  a  phenomenological 
sense.  In  particular,  we  would  prefer  not  to  require  a  large  number  of  additional 
parameters  to  achieve  a  reliable  representation  of  the  rheology  in  realistic  circumstances. 

A  carefully  designed  feedback  between  modeling  and  experimentation  appears  to  be  the 
appropriate  strategy  to  achieve  this  goal. 
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